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TECHNICAL NOTES 



Optimal short scan convolution reconstruction for fanbeam CT 

Dennis L Parker 

Department of Radiation Oncology, University of California at San Francisco, San Francisco, California 
94143 

(Received 22 June 1981; accepted for publication 23 October 1981) 

The problem of using a divergent fan beam convolution reconstruction algorithm in conjunction 
with a minimal complete (180° plus the fan angle) data set is reviewed. It is shown that by proper 
weighting of the initial data set, image quality essentially equivalent to the quality of 
reconstructions from 360? data sets is obtained. The constraints on the weights are that the sum of 
the two weights corresponding to the same line-integral must equal one, in regions of no data the 
weights must equal zero, and the weights themselves as well as the gradient of the weights must be 
continuous over the full 360*. After weighting the initial data with weights that satisfy these 
constraints, image reconstruction can be conveniently achieved by using the standard (hardwired 
if available) convolver and backprojector of the specific scanner. 



I. INTRODUCTION 

In the development of fast, special purpose CT scanners, a 
reduction in scanning time is also obtained by collecting only 




Fig. 1. Scanning geometry utilized in conventional fanbeam scanners. 



the minimum number of projections (views) required for im- 
age reconstruction. It is known, for example, that a scanner 
with a parallel beam geometry (such as a translate/rotate 
scanner) requires projections for 180 0 . 1 It has also been 
shown that data collected from fanbeam geometries can be 
rebinned or reorganized into parallel ray data sets for recon- 
struction with parallel geometry algorithms. 2 It is easy to 
show that a complete parallel ray equivalent data set can be 
obtained from a set of divergent ray projections taken over 
180° plus the angle of the divergent fan beam. It is logical to 
refer to this set of divergent projections which covers 180° 
plus the fan angle as the "minimal complete data set". It is 
the minimum set of equally spaced projection measurements 
which can be used in conventional convolution type recon- 
struction algorithms. In this note we show that with proper 
weighting of the input data set it is possible to reconstruct the 
minimal complete data set using conventional divergent 
convolution reconstruction algorithms (i.e., with no rebin- 
ning into parallel ray data sets). It is therefore possible to 
utilize conventional array processors and standard hard- 
wired backprojectors in performing the reconstruction of 
data from short scans. 

In a paper by Naparstek, 3 several methods of short scan 
fan beam reconstruction techniques were derived and as- 
sessed. However, those algorithms which generated images 
of reasonable quality also required more projections than a 
minimal complete set. The more successful algorithms uti- 
lized a smooth weighting of the extra projections in a manner 
similar to that used in standard CT overscanning for artifact 
reduction. 4 In this note, we show that by a simple modifica- 
tion of one of the techniques of Naparstek, good quality im- 
ages can be generated from a minimal complete short scan 
data set. 

II. ALGORITHM DEVELOPMENT 

Figure 1 illustrates the geometry of fan beam data collec- 
tion where 0 is the angle of rotation of the source to axis 
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Fig. 2. Outline of sampling region as a 
function of a (the horizontal coordinate) 
and 0 (the vertical coordinate). The upper 
triangle is geometrically the same (i.e., it 
samples the same set of line integrals) as 
the lower triangle. 



ANGLE WITHIN 
PROJ ECTION 



segment and a is the relative angular position of an individ- 
ual line integral from the source to a detector element If 
g(aj?) is the projection measurement (line integral) along the 
segment specified by (aJ3 ) and if p(t t # ) is the corresponding 
parallel ray projection measurement where (f ,t? ) are defined 
in the figure, it is evident that 

p{t$) = g[afi), (1) 
when 

# = 0 + a 
and 

t = D sin a. 
The symmetry over rotation requires that 

P(t,&)=p{-t 9 # + ir), (2) 
which is equivalent to 

g(a£) = g{ - a,ir +0 + 2a). (3) 
Thus, if data is collected in the rectangle, 
— 8<a <8, 

O<0<ir+28 9 (4) 

where 25 is the fan angle denned in Fig. 1, then the data in 
the triangle, 

0</3<28-2a, 

will be positionally redundant with the data in the region: 

ir-2a<0<ir + 26. 

This is illustrated in Fig. 2. The set of line integrals specified 
by the upper triangle is resampled in the lower triangle as 
required by Eq. (3). In contrast, the set of line integrals speci- 
fied by the parallelogram region bounded by the two trian- 
gles is only sampled once. 



Because only part of the region in the [aft) diagram is 
doubly scanned, it is expected that direct application of any 
direct divergent fanbeam reconstruction algorithm which it- 
self is based on total double scanning (360") would yield a 
distorted image. That this is true is illustrated using the sim- 
ulated data of Fig. 3. The mathematical phantom of Fig. 3(a) 
yields the divergent sinogram of Fig. 3(b) for a third genera- 
tion (rotate-rotate) scanner with I2l detectors at 1/3-degree 
angular increments for a fan angle of 40° and 360 views taken 
at one-degree projection increments. The data has been con- 
volved with a simple divergent algorithm kernel 5 : 



ir 2 sin 2 (£da)' 
q { =0,i,eveni>tO 



tea 2 



(5) 



When the entire (360°) data set is reconstructed using a 
standard backprojection algorithm, the image of Fig. 3(c) is 
obtained. If, however, the case of a minimal complete data 
set is imitated by considering only the first 220 projections in 
the standard backprojection algorithm, the image of Fig. 
3(d) is obtained. The variation in intensity over the image is 
surprisingly smooth and the image itself is very 
recognizable. 

It might be expected that the smooth variation could be 
removed by setting the data from one of the redundant trian- 
gles to zero. It was shown by Naparstek 3 that the streak 
artifacts which resulted from the zeroed boundary of the 
triangle were sufficient to completely mask the image. Of the 
other algorithms proposed by Naparstek, the most success- 
ful ones required additional data (beyond the minimal com- 
plete short scan data set). This additional data was then used 
to feather smoothly to zero the abrupt discontinuity at the 
parallelogram boundary. In addition to requiring data be- 
yond that of a minimal complete data set, the algorithms also 
generally required a more complicated convolution kernel 

The logical extension, which was not formally proposed 
or attempted by Naparstek, is to weight the data of the mini- 
mal complete data set in such a manner that the discontinu- 
ity is as uniformly distributed as possible. The constraints on 
the set of weights can be developed as follows: The weighted 
data set is defined as g*{aft ) and is given by 

flpfi^aiaftWaft), (6) 
where 

a>{afi) + o>(-a,ir + 0 + 2a)=\. (7) 
The minimization of the discontinuities at the various bor- 
ders of the triangles is accomplished by the requirements 

0(a,O) = 6>(a t ir + 23) = 0, (8) 

<o{a t 28 -2a) = &(a,ir — 2a) = l , (9) 



da(afi) 



0*=26-2a 



dP \0~ir-2a 

_ dco[a£) 



= 0. 



BP 
da(a&) \ 

da da 
One set of weights which satisfies all of the above con- 
straints is given by 



0*=> ir— 2a 



(10) 

(ii) 
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Fig. 3. (a) The mathematical 
phantom used in the simula- 
tions. See text for description, 
(b) The set of projection mea- 
surements in the same format 
as Fig. 2 after convolving with 
the reconstruction kernel. See 
text for geometry of simulated 
scanner. Because line integrals 
of negligible width were calcu- 
lated, the artifacts evident are 
due to aliasing, (c) The recon- 
struction of the data of Fig. 3(b) 
using a standard divergent 
backprojection algorithm and 
projection data over 360 *. (d) 
The reconstruction of the data 
of Fig. 3(b| using the same algo- 
rithm but with projection data 
over only 220 




a[a& ) = sin 2 ( i O<0<25 - 2a , 

V4 6 — a/ 

co[a£ ) = 1, 2<5 - 2a<fi<v - 2a , 

(o{afi ) = sin 2 [— ir + ^*~~ & \ ^ _ 2a<J3<ir + 26 . (12) 
\4 a / 



III. RESULTS 

In Fig. 4(a) we show the sinogram from Fig. 3(b) weighted 
by the functions of Eq. (12). Figure 4(b) is the reconstruction 
of the data from 4(a) using a standard divergent fanbeam 
backprojection algorithm. [This is the same algorithm as 
used in Figs. 3(c) and (d).] It is evident in comparing Figs. 3(c) 
and 4(b) that there is no appreciable difference in image qual- 
ity between the full 360° reconstruction and the 220°, weight- 
ed reconstruction. 



IV. CONCLUSION 

The technique which has been presented has demonstrat- 
ed the ability to reconstruct images from a minimal complete 
data set (180 9 plus the fan angle) with no evident loss in 
image quality over that obtained from reconstruction of the 
full 360 * scan data. The technique requires the smooth 
weighting of data in double-scanned regions while the 
weights and their derivatives are required to be continuous 
at the boundaries. After the weighting is performed, stan- 
dard divergent convolution/backprojection algorithms are 
used to reconstruct the image. It is therefore possible to 
maintain the hardware-related reconstruction speed of stan- 
dard convolvers and backprojectors while significantly re- 
ducing the required number of scanned projections, thereby 
decreasing the associated scan time. In this manner artifacts 
related to the duration of the scan, such as motion artifacts 
or electronic drifts in the data acquisition system, etc., can be 
minimized. 
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Fig. 4. (a) The projection data over 220 ' where the input data has been 
weighted by the weights given in the text and then convolved with the same 
algorithm as that of Fig. 3(b). (b) The reconstruction of the data of Fig. 4{a) 
using the same algorithm as Figs. 3(c) and (d|. 
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3 Algorithms for Reconstruction 
with Nondiffracting Sources 



In this chapter we will deal with the mathematical basis of tomography with 
nondiffracting sources. We will show how one can go about recovering the 
image of the cross section of an object from the projection data. In ideal 
situations, projections are a set of measurements of the integrated values of 
some parameter of the object— integrations being along straight lines through 
the object and being referred to as line integrals. We will show that the key to 
tomographic imaging is the Fourier Slice Theorem which relates the 
measured projection data to the two-dimensional Fourier transform of the 
object cross section. 

This chapter will start with the definition of line integrals and how they are 
combined to form projections of an object. By finding the Fourier transform 
of a projection taken along parallel lines, we will then derive the Fourier Slice 
Theorem. The reconstruction algorithm used depends on the type of 
projection data measured; we will discuss algorithms based on parallel beam 
projection data and two types of fan beam data. 

3.1 Line Integrals and Projections 

A line integral, as the name implies, represents the integral of some 
parameter of the object along a line. In this chapter we will not concern 
ourselves with the physical phenomena that generate line integrals, but a 
typical example is the attenuation of x-rays as they propagate through 
biological tissue. In this case the object is modeled as a two-dimensional (or 
three-dimensional) distribution of the x-ray attenuation constant and a line 
integral represents the total attenuation suffered by a beam of x-rays as it 
travels in a straight line through the object. More details of this process and 
other examples will be presented in Chapter 4. 

We will use the coordinate system defined in Fig. 3.1 to describe line 
integrals and projections. In this example the object is represented by a two- 
dimensional function /(x, y) and each line integral by the (0, /) parameters. 

The equation of line AB in Fig. 3. 1 is 

x cos 0+y sin 0 = / (1) 




Fig. 3.1: An object, f(x, y)> 
and its projection, P 8 (t x ), are 
shown for an angle of 0. (From 
lKak79l) 



and we will use this relationship to define line integral P e (t) as 

Pe(*)=\ f(x 9 y)ds. 

J (M line 

Using a delta function, this can be rewritten as 

Pe(t)= P (" /(*> y)&(* cos 0+>> sin d-t) dx dy. 

J _ 0) J — CO 

The function P 9 (t) is known as the Radon transform of the function /(x, 
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Fig. 3.2: Parallel projections are A projection is formed by combining a set of line integrals. The simplest 

taken by measuring a set of projection is a collection of parallel ray integrals as is given by P e (0 for a 

parallel rays for °™™*"°{*,, constant 9. This is known as a parallel projection and is shown in Fig. 3.2. It 
different angles. (From [Ros82J.) , , - j j ^ * 

could be measured, for example, by moving an x-ray source and detector 

along parallel lines on opposite sides of an object. 

Another type of projection is possible if a single source is placed in a fixed 
position relative to a line of detectors. This is shown in Fig. 3.3 and is known 
as a fan beam projection because the line integrals are measured along fans. 

Most of the computer simulation results in this chapter will be shown for 
the image in Fig. 3.4. This is the well-known Shepp and Logan [She74] 
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P(t,B 2 ) 




Fig. 3.3: A fan beam projection 
is collected if all the rays meet in 
one location. (From [Ros82]J 



"head phantom," so called because of its use in testing the accuracy of 
reconstruction algorithms for their ability to reconstruct cross sections of the 
human head with x-ray tomography. (The human head is believed to place the 
greatest demands on the numerical accuracy and the freedom from artifacts of 
a reconstruction method.) The image in Fig. 3.4(a) is composed of 10 
ellipses, as illustrated in Fig. 3.4(b). The parameters of these ellipses are 
given in Table 3.1. 

A major advantage of using an image like that in Fig. 3.4(a) for computer 
simulation is that now one can write analytical expressions for the 
projections. Note that the projection of an image composed of a number of 
ellipses is simply the sum of the projections for each of the ellipses; this 
follows from the linearity of the Radon transform. We will now present 
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-1.0-0.8-0.6-0.^ -0.2 0.0 0.2 0.*« 0.6 0.8 1.0 

x 

(b) 

Fig. 3,4: The Shepp and Logan expressions for the projections of a single ellipse. Let/(x, y) be as shown in 
"head phantom" is shown in (a). pjg t 3.5(a), i.e., 
Most of the computer simulated 9 ' 

results in this chapter were 
generated using this phantom. 
The phantom is a superposition 

of 10 ellipses, each with a size /(•*» >0 

and magnitude as shown in (b). 
(From [Ros82].) 



p for — + — < 1 (inside the ellipse) 
A 2 Br 

< 0 otherwise (outside the ellipse) . 
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Fig. 3.5: (a) An analytic 
expression is shown for the 
projection of an ellipse. For 
computer simulations a projection 
can be generated by simply 
summing the projection of each 
individual ellipse, (b) Shown here 
is an ellipse with its center located 
at (JC|, y\) and its major axis 
rotated by a. (From fRos82J.) 




X 
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Fig. 3.5: Continued. 
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Table 3.1: Summary of parameters for tomography simulations. 



Center 


Major 


Minor 


Rotation 


Refractive 


Coordinate 


Axis 


Axis 


Angle 


Index 


(0, 0) 


0.92 


0.69 


90 


2.0 


(0, -0.0184) 


0.874 


0.6624 


90 


-0.98 


(0.22, 0) 


0.31 


0.11 


72 


-0.02 


(-0.22, 0) 


0.41 


0.16 


108 


-0.02 


(0, 0.35) 


0.25 


0.21 


90 


0.01 


(0, 0.1) 


0.046 


0.046 


0 


0.01 


(0, -0.1) 


0.046 


0.046 


0 


0.01 


(-0.08, -0.605) 


0.046 


0.023 


0 


0.01 


(0, -0.605) 


0.023 


0.023 


0 


0.01 


(0.06, -0.605) 


0.046 


0.023 


90 


0.01 
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It is easy to show that the projections of such a function are given by 

f^V^W^ for \t\*a(d) 

(O \t\>a(6) 

where a\0) = A 2 cos 2 0 + B z sin 2 0. Note that a(6) is equal to the projection 
half-width as shown in Fig. 3.5(a). 

Now consider the ellipse described above centered at (x u y\) and rotated 
by an angle a as shown in Fig. 3.5(b). Let P'(0, t) be the resulting 
projections. They are related to P B {t) in (5) by 

P e (t) = P e -a(t-s cos (7 - 0)) (6) 

where s = \ix] + y] and y = tan" 1 (j\/xi). 



3.2 The Fourier Slice Theorem 



We derive the Fourier Slice Theorem by taking the one-dimensional 
Fourier transform of a parallel projection and noting that it is equal to a slice 
of the two-dimensional Fourier transform of the original object. It follows 
that given the projection data, it should then be possible to estimate the object 
by simply performing a two-dimensional inverse Fourier transform. 

We start by defining the two-dimensional Fourier transform of the object 
function as 

F(u,v)=\ m T f(x y y)e-^ ux ^dxdy. (7) 

Likewise define a projection at an angle 0, P 0 (t), and its Fourier transform by 

S,(>v)=r P e (t)e-J 2 ™' dt. (8) 

The simplest example of the Fourier Slice Theorem is given for a 
projection at 0 = 0. First, consider the Fourier transform of the object along 
the line in the frequency domain given by v = 0. The Fourier transform 
integral now simplifies to 

F(u>0)=[" f f(x y y)e^ 2rWC dx dy (9) 

but because the phase factor is no longer dependent on y we can split the 
integral into two parts, 

f( Wj o)= J 00 [ ry(*> y) d y\ e ' J2rux dx - < 10 > 

From the definition of a parallel projection, the reader will recognize the term 
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in brackets as the equation for a projection along lines of constant x or 

P,.o(x)=t" f(x,y)dy. (H) 
Substituting this in (10) we find 

F(«, 0)=f Pt-oMe-' 1 "" dx. (12) 

J _ 00 

The right-hand side of this equation represents the one-dimensionai Fourier 
transform of the projection P e=0 \ thus we have the following relationship 
between the vertical projection and the 2-D transform of the object function: 

F(u,0) = S^o(«). 03) 

This is the simplest form of the Fourier Slice Theorem. Clearly this result 
is independent of the orientation between the object and the coordinate 
system. If, for example, as shown in Fig. 3.6 the (/, s) coordinate system is 
rotated by an angle 0, the Fourier transform of the projection defined in (1 1) 
is equal to the two-dimensional Fourier transform of the object along a line 
rotated by 6. This leads to the Fourier Slice Theorem which is stated as 
[Kak85]: 

The Fourier transform of a parallel projection of an image f(x, y) 
taken at angle $ gives a slice of the two-dimensional transform, 
F(w, v), subtending an angle 0 with the w-axis. In other words, 
the Fourier transform of P 9 (t) gives the values of F{u> v) along 
line BB in Fig. 3.6. 



Fig. 3.6: The Fourier Siice 
Theorem relates the Fourier 
transform of a projection to the 
Fourier transform of the object 
along a radial line. (From 
[Pan83].) 




ALGORITHMS FOR RECONSTRUCTION WITH NONDIFFRACTING SOURCES 57 



The derivation of the Fourier Slice Theorem can be placed on a more solid 
foundation by considering the (/, s) coordinate system to be a rotated version 
of the original (x, y) system as expressed by 

[/I f cos 0 sin 0 1 [xl /tA . 
,J = [-sin0cos*JUJ * (H) 

In the (t, s) coordinate system a projection along lines of constant t is written 

Pe(t)=\~ f«,s)ds (15) 

J - CO 

and from (8) its Fourier transform is given by 

S 9 (w)=T P t {f)e-* u « dt. (8) 

J - GO 

Substituting the definition of a projection into the above equation we find 
S e (w) = J" 5) <fc| £T' 2 ™" dL (16) 

This result can be transformed into the (x> y) coordinate system by using the 
relationships in (14), the result being 

5 fl (w)=r P f(x, y)e-j 2 ™( x ™ $ +y™v dxdy. (17) 

v — CO J — OO 

The right-hand side of this equation now represents the two-dimensional 
Fourier transform at a spatial frequency of (u = w cos 0, v = w sin 0) or 

S e (w)=F(w, 0)=F(w cos 0, w sin 0). (18) 

This equation is the essence of straight ray tomography and proves the 
Fourier Slice Theorem. 

The above result indicates that by taking the projections of an object 
function at angles 0 u 0i, • * • , 0* and Fourier transforming each of these, we 
can determine the values of F(u, v) on radial lines as shown in Fig. 3.6. If an 
infinite number of projections are taken, then F(u, v) would be known at all 
points in the wi/-plane. Knowing F(w, v), the object function /(x, y) can be 
recovered by using the inverse Fourier transform: 

/(x, y)= P P F(u, i/)e>2*<"*+»» du dv. (19) 

V — CD J - CO 

If the function /(x, y) is bounded by -A/2 < x < A/2 and - A/2 < y < 
A/2, for the purpose of computation (19) can be written as 

fiX, y) = ^- 2 2 £ F efriimtAWntAW (20) 
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for 



Fig. 3.7: Collecting projections 
of the object at a number of 
angles gives estimates of the 
Fourier transform of the object 
along radial lines. Since an FFT 
algorithm is used for 
transforming the data, the dots 
represent the actual location of 
estimates of the object's Fourier 
transform. (From fPan83JJ 



A A A A 

~2 <X< 2 ^ "2 <y< 2 



(21) 



Since in practice only a finite number of Fourier components will be known, 
we can write 



1 



N/l 



N/2 



m= -N/l n=-N/2 



2*((m/A)x+(n/A)y) 



for 



A A A A 

— <x<— and <y<— 

2 2 2 2 



(22) 



(23) 



where we arbitrarily assume N to be an even integer. It is clear that the spatial 
resolution in the reconstructed picture is determined by N. Equation (22) can 
be rapidly implemented by using the fast Fourier transform (FFT) algorithm 
provided the N 2 Fourier coefficients F(m/A, n/A) are known. 

In practice only a finite number of projections of an object can be taken. In 
that case it is clear that the function F(u, v) is only known along a finite 
number of radial lines such as in Fig. 3.7. In order to be able to use (22) one 
must then interpolate from these radial points to the points on a square grid. 
Theoretically, one can exactly determine the N 2 coefficients required in (22) 
provided as many values of the function F(u, v) are known on some radial 
lines [Cro70]. This calculation involves solving a large set of simultaneous 
equations often leading to unstable solutions. It is more common to determine 
the values on the square grid by some kind of nearest neighbor or linear 
interpolation from the radial points. Since the density of the radial points 
becomes sparser as one gets farther away from the center, the interpolation 
error also becomes larger. This implies that there is greater error in the 




frcqueney domain 
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calculation of the high frequency components in an image than in the low 
frequency ones, which results in some image degradation. 

3.3 Reconstruction Algorithms for Parallel Projections 

The Fourier Slice Theorem relates the Fourier transform of a projection to 
the Fourier transform of the object along a single radial. Thus given the 
Fourier transform of a projection at enough angles the projections could be 
assembled into a complete estimate of the two-dimensional transform and 
then simply inverted to arrive at an estimate of the object. While this provides 
a simple conceptual model of tomography, practical implementations require 
a different approach. 

The algorithm that is currently being used in almost all applications of 
straight ray tomography is the filtered backprojection algorithm. It has been 
shown to be extremely accurate and amenable to fast implementation and will 
be derived by using the Fourier Slice Theorem. This theorem is brought into 
play by rewriting the inverse Fourier transform in polar coordinates and 
rearranging the limits of the integration therein. The derivation of this 
algorithm is perhaps one of the most illustrative examples of how we can 
obtain a radically different computer implementation by simply rewrit- 
ing the fundamental expressions for the underlying theory. 

In this chapter, derivations and implementation details will be presented for 
the backprojection algorithms for three types of scanning geometries, parallel 
beam, equiangular fan beam, and equispaced fan beam. The computer 
implementation of these algorithms requires the projection data to be sampled 
and then filtered. Using FFT algorithms we will show algorithms for fast 
computer implementation. Before launching into the mathematical deriva- 
tions of the algorithms, we will first provide a bit of intuitive rationale behind 
the filtered backprojection type of approach. If the reader finds this 
presentation excessively wordy, he or she may go directly to Section 3.3.2. 

3.3.1 The Idea 

The filtered backprojection algorithm can be given a rather straightforward 
intuitive rationale because each projection represents a nearly independent 
measurement of the object. This isn't obvious in the space domain but if the 
Fourier transform is found of the projection at each angle then it follows 
easily by the Fourier Slice Theorem. We say that the projections are nearly 
independent (in a loose intuitive sense) because the only common information 
in the Fourier transforms of the two projections at different angles is the dc 
term. 

To develop the idea behind the filtered backprojection algorithm, we note 
that because of the Fourier Slice Theorem the act of measuring a projection 
can be seen as performing a two-dimensional filtering operation. Consider a 
single projection and its Fourier transform. By the Fourier Slice Theorem, 
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Fig. 3.8: This figure shows the 
frequency domain data available 
from one projection, (a) is the 
ideal situation. A reconstruction 
could be formed by simply 
summing the reconstruction from 
each angle until the entire 
frequency domain is filled. What 
is actually measured is shown in 
(b). As predicted by the Fourier 
Slice Theorem, a projection gives 
information about the Fourier 
transform of the object along a 
single line. The filtered 
backprojection algorithm takes 
the data in (b) and applies a 
weighting in the frequency 
domain so that the data in (c) are 
an approximation to those in (a). 



this projection gives the values of the object's two-dimensional Fourier 
transform along a single line. If the values of the Fourier transform of this 
projection are inserted into their proper place in the object's two-dimensional 
Fourier domain then a simple (albeit very distorted) reconstruction can be 
formed by assuming the other projections to be zero and finding the two- 
dimensional inverse Fourier transform. The point of this exercise is to show 
that the reconstruction so formed is equivalent to the original object's Fourier 
transform multiplied by the simple filter shown in Fig. 3.8(b). 

What we really want from a simple reconstruction procedure is the sum of 
projections of the object filtered by pie-shaped wedges as shown in Fig. 
3.8(a). It is important to remember that this summation can be done in either 
the Fourier domain or in the space domain because of the linearity of the 
Fourier transform. As will be seen later, when the summation is carried out in 
the space domain, this constitutes the backprojection process. 

As the name implies, there are two steps to the filtered backprojection 
algorithm: the filtering part, which can be visualized as a simple weighting of 
each projection in the frequency domain, and the backprojection part, which 
is equivalent to finding the elemental reconstructions corresponding to each 
wedge filter mentioned above. 

The first step mentioned above accomplishes the following: A simple 
weighting in the frequency domain is used to take each projection and 
estimate a pie-shaped wedge of the object's Fourier transform. Perhaps the 
simplest way to do this is to take the value of the Fourier transform of the 
projection, S e (w), and multiply it by the width of the wedge at that frequency. 
Thus if there are K projections over 180° then at a given frequency w, each 
wedge has a width of 2-k\w\/K. Later when we derive the theory more 
rigorously, we will see that this factor of |w| represents the Jacobian for a 
change of variable between polar coordinates and the rectangular coordinates 
needed for the inverse Fourier transform. 

The effect of this weighting by 2r\w\/K is shown in Fig. 3.8(c). 
Comparing this to that shown in (a) we see that at each spatial frequency, u>, 
the weighted projection, (2tt| w\/K)S e (w), has the same "mass" as the pie- 
shaped wedge. Thus the weighted projections represent an approximation to 
the pie-shaped wedge but the error can be made as small as desired by using 
enough projections. 

The final reconstruction is found by adding together the two-dimensional 
inverse Fourier transform of each weighted projection. Because each 
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projection only gives the values of the Fourier transform along a single line, 
this inversion can be performed very fast. This step is commonly called a 
backprojection since, as we will show in the next section, it can be perceived 
as the smearing of each filtered projection over the image plane. 
The complete filtered backprojection algorithm can therefore be written as: 

Sum for each of the K angles, 0, between 0 and 180° 
Measure the projection, P 6 (t) 
Fourier transform it to find S e (w) 
Multiply it by the weighting function 2ir\w\/K 
Sum over the image plane the inverse Fourier transforms of the 
filtered projections (the backprojection process). 

There are two advantages to the filtered backprojection algorithm over a 
frequency domain interpolation scheme. Most importantly, the reconstruction 
procedure can be started as soon as the first projection has been measured. 
This can speed up the reconstruction procedure and reduce the amount of data 
that must be stored at any one time. To appreciate the second advantage, the 
reader must note (this will become clearer in the next subsection) that in the 
filtered backprojection algorithm, when we compute the contribution of each 
filtered projection to an image point, interpolation is often necessary; it turns 
out that it is usually more accurate to carry out interpolation in the space 
domain, as part of the backprojection or smearing process, than in the 
frequency domain. Simple linear interpolation is often adequate for the 
backprojection algorithm while more complicated approaches are needed for 
direct Fourier domain interpolation [Sta81]* 

In Fig. 3.9(a) we show the projection of an ellipse as calculated by (5). To 
perform a reconstruction it is necessary to filter the projection and then 
backproject the result as shown in Fig. 3.9(b). The result due to backproject- 
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Fig, 3,9: A projection of an 
ellipse is shown in (a), (b) shows 
the projection after it has been 
filtered in preparation for 
backprojection. 



Fig. 3.10: The result of 
backprojecting the projection in 
Fig, 3.9 is shown here, (a) shows 
the result of backprojecting for a 
single angle, (b) shows the effect 
of backprojecting over 4 angles, 
(c) shows 64 angles, and (d) 
shows 512 angles. 



ing one projection is shown in Fig. 3.10. It takes many projections to 
accurately reconstruct an object; Fig. 3.10 shows the result of reconstructing 
an object with up to 512 projections. 

3.3.2 Theory 

We will first present the backprojection algorithm for parallel beam 
projections. Recalling the formula for the inverse Fourier transform, the 
object function, f(x, y), can be expressed as 



/(a:, y)= T f" F(u 9 v )e^ ux ^y ) du dv. 

J _ oe « — 00 



(24) 



Exchanging the rectangular coordinate system in the frequency domain, (w, 
u), for a polar coordinate system, (w, 0), by making the substitutions 



u - w cos 9 
v = w sin 6 

and then changing the differentials by using 

du dv = w dw d6 



(25) 
(26) 

(27) 



(a*. t^pysdJon 
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we can write the inverse Fourier transform of a polar function as 

fix, y) = f T (°° F(w, 6)e> 2vw(x «* e+y e) w dw dd. (28) 
Jo Jo 

This integral can be split into two by considering 8 from 0° to 180° and then 
from 180° to 360°, 

Ax, y)=\* P F(w, fl^^cxcos 9 j w m 
Jo Jo 

+ (' [" F(w, e+m°)e> Uw[xcos < i+me)+} ' sin(9+, * 0 '' )1 w dw dd, 

Jo *0 



(29) 



and then using the property 

F(w, 0+18O o )=F(-w, 0) (30) 
the above expression for f(x, y) may be written as 

/( X , y)= j* |^ j°° F(w, 0)\w\eJ Uwt c/wj ctf. (31) 

Here we have simplified the expression by setting 

/=jccos 0+y sin 0. (32) 

If we substitute the Fourier transform of the projection at angle 0, S e (w), for 
the two-dimensional Fourier transform F(w t 0), we get 

f(x, y)= j* |^ j"^ S B (w)\w\e^ </>vj d0. (33) 

This integral in (33) may be expressed as 

f( x , y) = (' Q,(x cos 0 + y sin 0) <*0 (34) 

J o 

where 

Qe(t)=[° S d (w)\w\e' Uwt dw. (35) 

This estimate of /(*, >>), given the projection data transform S*(w), has a 
simple form. Equation (35) represents a filtering operation, where the 
frequency response of the filter is given by | w\; therefore Q $ (w) is called a 
"filtered projection." The resulting projections for different angles 0 are then 
added to form the estimate of /(x, y). 

Equation (34) calls for each filtered projection, Q $ , to be "backpro- 
jected. " This can be explained as follows. To every point (x, y) in the image 
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plane there corresponds a value of / = x cos 6 + y sin 6 for a given value of 
0, and the filtered projection Q e contributes to the reconstruction its value at / 
( = x cos 6 + y sin 0). This is further illustrated in Fig. 3.11. It is easily 
shown that for the indicated angle d, the value of t is the same for all (x 9 y) on 
the line LM. Therefore, the filtered projection, Q e , will make the same 
contribution to the reconstruction at all of these points. Therefore, one 
could say that in the reconstruction process each filtered projection, Q 0y is 
smeared back, or backprojected, over the image plane. 

The parameter w has the dimension of spatial frequency. The integration in 
(35) must, in principle, be carried out over all spatial frequencies. In practice 
the energy contained in the Fourier transform components above a certain 
frequency is negligible, so for all practical purposes the projections may be 
considered to be bandlimited. If W is a frequency higher than the highest 
frequency component in each projection, then by the sampling theorem the 
projections can be sampled at intervals of 



Fig. 3.11: Reconstructions are 
often done using a procedure 
known as backprojection. Here a 
filtered projection is smeared 
back over the reconstruction 
plane along lines of constant t. 
The filtered projection at a point 
t makes the same contribution to 
all pixels along the line LM in the 
x-y plane. (From fRos82J.) 



T=- 



2W 



(36) 



without introducing any error. If we also assume that the projection data are 
equal to zero for large values of \t\ then a projection can be represented as 



Pe(mT), m = - 



-N 



N 
■>2~ l 



(37) 



for some (large) value of TV. An FFT algorithm can then be used to 




approximate the Fourier transform S e (w) of a projection by 

«„>.» (. 2 -g)-±x (^) ™ 

Given the samples of a projection, (38) gives the samples of its Fourier 
transform. The next step is to evaluate the "modified projection" Q e (t) 
digitally. Since the Fourier transforms S e (w) have been assumed to be 
bandlimited, (35) can be approximated by 

<MO=r S $ {w)\w\eJ Uw dw (39) 
J - w 



2W / 2^\| 2W 

77 S s »{ m l^)\ m 17 



e jlrm{2W/N)t 



provided TV is large enough. Again, if we want to determine the projections 
Q e (t) for only those / at which the projections P 0 (t) are sampled, we get 



/ k \ (2W\ aw o / 2W\ 



2W 
/n — 
N 



e jlx{mk/N) (41) 

m = - /v/2 N ' 

k=-N/2 9 -1,0, 1, •••,N/2. (42) 

By the above equation the function (?*(/) at the sampling points of the 
projection functions is given (approximately) by the inverse DFT of the 
product of S e (m(2 W/N)) and \m(2 W/N)\ . From the standpoint of noise in 
the reconstructed image, superior results are usually obtained if one 
multiplies the filtered projection, S e (2W/N)\m(2W/N)\, by a function 
such as a Hamming window [Ham77]: 

(£M5)JL*(-5) 

H^^^je^W'W (43) 

where H(m(2W/N)) represents the window function used. The purpose of 
the window function is to deemphasize high frequencies which in many cases 
represent mostly observation noise. By the familiar convolution theorem for 
the case of discrete transforms, (43) can be written as 

<"> 

where * denotes circular (periodic) convolution and where 0(A:/2 W) is the 
inverse DFT of the discrete function \m(2W/N)\H(m(21V/N)), m = 
-N/2, -1, 0, 1, '--,N/2. 



n 

I 2W 
• m — 
N 
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Clearly at the sampling points of a projection, the function Q e (t) may be 
obtained either in the Fourier domain by using (40), or in the space domain by 
using (44). The reconstructed picture /(x, y) may then be obtained by the 
discrete approximation to the integral in (34), i.e., 

/(*» y) = p 2 W* cos ft +y sin w (45) 

where the A' angles 0, are those for which the projections P $ (t) are known. 

Note that the value of x cos 0, + 7 sin 0/ in (45) may not correspond to one 
of the values of t for which Q 9i is determined in (43) or in (44). However, Q 9i 
for such t may be approximated by suitable interpolation; often linear 
interpolation is adequate. 

Before concluding this subsection we would like to make two comments 
about the filtering operation in (35). First, note that (35) may be expressed in 
the /-domain as 

Q*{t)=\Pt{*)p(f-a)d* (46) 

where p(t) is nominally the inverse Fourier transform of the | w\ function in 
the frequency domain. Since |w| is not a square integrable function, its 
inverse transform doesn't exist in an ordinary sense. However, one may 
examine the inverse Fourier transform of 

|w|«-M (47) 

as € -+ 0. The inverse Fourier transform of this function, denoted by p&\ is 
given by 

This function is sketched in Fig. 3.12. Note that for large t we get p € (t) 
- 1/(2tt/) 2 . 

Now our second comment about the filtered projection in (35): This 
equation may also be written as 

&(/)= j2irwS 0 (w) ^ sgn (w)J e' u " dw (49) 



where 



, . f 1 for w>0 /cm 
^ forw<0 . ( 50 > 



By the standard convolution theorem, this equation may be expressed as 

Q e (t) = {IFT of j2irwS e (w)} * {IFT of sgn (w)} (51) 
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Fig. 3.12: An approximation to 
the impulse response of the ideal 
backprojection filter is shown 
here, (From [Ros82].) 



where the symbol * denotes convolution and the abbreviation IFT stands for 
inverse fast Fourier transform. The IFT of j2irwSe(w) is (d/dt)P $ (t) while 
the IFT of ( - j/2w) sgn (w) is \/t. Therefore, the above result may be written 
as 



1 dP 9 {t) 



(52) 



= Hilbert Transform of 



3P«(0 
dt 



(53) 



where, expressed as a filtering operation, the Hilbert Transform is usually 
defined as the following frequency response: 



«<»>-[-*: :<o. 



(54) 
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3.3.3 Computer Implementation of the Algorithm 

Let's assume that the projection data are sampled with a sampling interval 
of t cm. If there is no aliasing, this implies that in the transform domain the 
projections don't contain any energy outside the frequency interval (- W, 
W) where 

W=— cycles/cm. (55) 
2t 

Let the sampled projections be represented by Pq(Ict) where k takes integer 
values. The theory presented in the preceding subsection says that for each 
sampled projection Peikr) we must generate a filtered Q e (kr) by using the 
periodic (circular) convolution given by (40). Equation (40) is very attractive 
since it directly conforms to the definition of the DFT and, if N is 
decomposable, possesses a fast FFT implementation. However, note that (40) 
is only valid when the projections are of finite bandwidth and finite order. 
Since these two assumptions (taken together) are never strictly satisfied, 
computer processing based on (40) usually leads to interperiod interference 
artifacts created when an aperiodic convolution (required by (35)) is 
implemented as a periodic convolution. This is illustrated in Fig. 3.13. Fig. 
3.13(a) shows a reconstruction of the Shepp and Logan head phantom from 
1 10 projections and 127 rays in each projection using (40) and (45). Equation 
(40) was implemented with a base 2 FFT algorithm using 128 points. Fig. 
3.13(b) shows the reconstructed values on the horizontal line for y = 
-0.605. For comparison we have also shown the values on this line in the 
original object function. 

The comparison illustrated in Fig. 3.13(b) shows that reconstruction based 
on (42) and (45) introduces a slight "dishing" and a dc shift in the image. 
These artifacts are partly caused by the periodic convolution implied by (40) 
and partly by the fact that the implementations in (40) "zero out" all the 
information in the continuous frequency domain in the cell represented by m 
= 0, whereas the theory (eq. (35)) calls for such "zeroing out" to occur at 
only one frequency, viz. w = 0. The contribution to these artifacts by the 
interperiod interference can be eliminated by adequately zero-padding the 
projection data before using the implementations in (42) or (43). 

Zero-padding of the projections also reduces, but never completely 
eliminates, the contribution to the artifacts by the zeroing out of the 
information in the m = 0 cell in (40). This is because zero-padding in the 
space domain causes the cell size to get smaller in the frequency domain. (If 
Nfft points are used for performing the discrete Fourier transform, the size 
of each sampling cell in the frequency domain is equal to l/yV F FT r ) To 
illustrate the effect of zero-padding, the 127 rays in each projection in the 
preceding example were padded with 129 zeros to make the data string 256 
elements long. These data were transformed by an FFT algorithm and filtered 
with a |h>| function as before. The y = -0.605 line through the 
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Fig. 3.13: (a) This 
reconstruction of the Shepp and 
Logan phantom shows the 
artifacts caused when the 
projection data are not 
adequately zero-padded and FFTs 
are used to perform the filtering 
operation in the filtered 
backprojection algorithm. The 
dark regions at the top and the 
bottom of the reconstruction are 
the most visible artifacts here. 
This 128 x 128 reconstruction 
was made from 110 projections 
with 127 rays in each projection, 
(b) A numerical comparison of 
the true and the reconstructed 
values on the y - - 0.605 line. 
(For the location of this line see 
Fig. 3.4.) The "dishing" and the 
dc shift artifacts are quite evident 
in this comparison, (c) Shown 
here are the reconstructed values 
obtained on they = -0.605 line 
if the 127 rays in each projection 
are zero-padded to 2S6 points 
before using the FFTs. The 
dishing caused by interperiod 
interference has disappeared; 
however, the dc shift still 
remains. (From fRos82J.) 
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Fig. 3.13: Continued. 
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reconstruction is shown in Fig. 3.13(c), demonstrating that the dishing 
distortion is now less severe. 

We will now show that the artifacts mentioned above can be eliminated by 
the following alternative implementation of (35) which doesn't require the 
approximation used in the discrete representation of (40). When the highest 
frequency in the projections is finite (as given by (55)), (35) may be 
expressed as 



where 



where, again, 



Qe(t)= P S e (w)H(w)e^ wt dw 

J -00 

//(w)HHMw) 
b "W={o oth! 



\<w 

otherwise. 



(56) 
(57) 
(58) 



H(w), shown in Fig. 3.14, represents the transfer function of a filter with 
which the projections must be processed. The impulse response, h{t), of this 
filter is given by the inverse Fourier transform of H(w) and is 



/i(r)= P H (w)e + j Uwt dw 

J _ oo 

1 sin2*7/2T 1 / 



2t 2 2*7/27 



sin vt/2r\ 2 
■Kt/lr ) 



(59) 
(60) 
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Frequency (w) 



Fig. 3.14: The ideal filter 
response for the filtered 
backprojection algorithm is 
shown here. It has been 
bandlimited to l/2r. (From 
[Ros82U 



where we have used (55). Since the projection data are measured with a 
sampling interval of r, for digital processing the impulse response need only 
be known with the same sampling interval. The samples, h(nr), of h(t) are 
given by 

'l/4r 2 , n = 0 

0, n even (61) 

1 



h(nr) = 



n odd. 



This function is shown in Fig. 3.15. 

Since both P e (t) and h(t) are now bandlimited functions, they may be 
expressed as 



Pa(0= 2 Peikr) 

*= -a. 



sin IvWjt-kr) 
2vW(t-kr) 

sin 2vW{t-kT) 



(62) 



(63) 



*= -CO 



2vW{t-kr) 

By the convolution theorem the filtered projection (56) can be written as 

Q 6 (0=r PM'Mt-ndt'- («) 

J - GO 
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Fig. 3.15: 77* iwte rcpoiisr Substituting (62) and (63) in (64) we get the following result for the values of 
of the filter shown in Fig. 3.14 is the filtered projection at the sampling points: 
shown here. (From [Ros82JJ 

Qo(nT) = T J h(nr-kr)P e (kr). (65) 

In practice each projection is of only finite extent. Suppose that each P${kr) is 
zero outside the index range k = 0, • • N - 1. We may now write the 
following two equivalent forms of (65): 

Qe(nr) = T h(nr-kr)P B (kr\ n = 0> 1,2, ■ • N- 1 (66) 

or 



Q e (nT) = t ^ h(kr)P$(nr-kr) 9 * = 0, 1, 2, • ■•, N- 1. (67) 

* = -(N-l) 
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Fig. 3.16: TheDFTofthe 
bandlimited filter (broken line) 
and that of the ideal filter (solid 
line) are shown here. Notice the 
primary difference is in the dc 
component. (From [Ros82J.) 



These equations imply that in order to determine Qeinr) the length of the 
sequence h(nr) used should be from / = - (N - 1) to / = (N - 1). It is 
important to realize that the results obtained by using (66) or (67) aren't 
identical to those obtained by using (42). This is because the discrete Fourier 
transform of the sequence h (/it) with n taking values in a finite range [such as 
when n ranges from - (N - 1) to (N - 1)] is not the sequence \k[(2 W)/ 
N]|. While the latter sequence is zero at k = 0, the DFT of h{nr) with n 
ranging from -(N - 1) to (N - 1) is nonzero at this point. This is 
illustrated in Fig. 3.16. 

The discrete convolution in (66) or (67) may be implemented directly on a 
general purpose computer. However, it is much faster to implement it in the 
frequency domain using FFT algorithms. [By using specially designed 
hardware, direct implementation of (66) can be made as fast or faster than the 
frequency domain implementation.] For the frequency domain implementa- 
tion one has to keep in mind the fact that one can now only perform periodic 
(or circular) convolutions, while the convolution required in (66) is aperiodic. 
To eliminate the interperiod interference artifacts inherent to periodic 
convolution, we pad the projection data with a sufficient number of zeros. It 
can easily be shown [Jak76] that if we pad P B (kr) with zeros so that it is (27V 
- 1) elements long, we avoid interperiod interference over the N samples of 
Qe(kr). Of course, if one wants to use the base 2 FFT algorithm, which is 
most often the case, the sequences P e {kr) and h(kr) have to be zero-padded 
so that each is (2N - 1) 2 elements long, where (2N - 1) 2 is the smallest 
integer that is a power of 2 and that is greater than 2N - 1 . Therefore, the 
frequency domain implementation may be expressed as 



£M"T) = rxIFFT {[FFT P 9 {m) with ZP]x[FFT h(nr) with ZP]}, (68) 
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where FFT and IFFT denote, respectively, fast Fourier transform and inverse 
fast Fourier transform; ZP stands for zero-padding. One usually obtains 
superior reconstructions when some smoothing is also incorporated in (68). 
Smoothing may be implemented by multiplying the product of the two FFTs 
by a Hamming window. When such a window is incorporated, (68) may be 
rewritten as 

Q e (nr) = TXlFFT {[FFT P e (nr) with ZP] 

X [FFT h (nr) with ZP] X smoothing - window } . (69) 

After the filtered projections Qe(m) are calculated with the alternative 
method presented here, the rest of the implementation for reconstructing the 
image is the same as in the preceding subsection. That is, we use (45) for 
backprojections and their summation. Again for a given (x, y) and 0, the 
argument x cos 0, + y sin 0,- may not correspond to one of the kr at which Q Bi 
is known. This will call for interpolation and often linear interpolation is 
adequate. Sometimes, in order to eliminate the computations required for 
interpolation, preinterpolation of the functions Q$(t) is also used. In this 
technique, which can be combined with the computation in (69), prior to 
backprojection, the function Q $ (t) is preinterpolated onto 10 to 1000 times 
the number of points in the projection data. From this dense set of points one 
simply retains the nearest neighbor to obtain the value of Q $i at x cos 0/ + y 
sin 0/. A variety of techniques are available for preinterpolation [Sch73]. 

One method of preinterpolation, which combines it with the operations in 
(69), consists of the following: In (69), prior to performing the IFFT, the 
frequency domain function is padded with a large number of zeros. The 
inverse transform of this sequency yields the preinterpolated Q 9 . It was 
recently shown [Kea78] that if the data sequence contains "fractional" 
frequencies this approach may lead to large errors especially near the 
beginning and the end of the data sequence. Note that with preinterpolation 
and with appropriate programming, the backprojection for parallel projection 
data can be accomplished with virtually no multiplications. 

Using the implementation in (68), Fig. 3.17(b) shows the reconstructed 
values on the line y - - 0,605 for the Shepp and Logan head phantom. 
Comparing with Fig. 3.13(b), we see that the dc shift and the dishing have 
been eliminated. Fig. 3.17(a) shows the complete reconstruction. The 
number of rays used in each projection was 127 and the number of projections 
100. To make convolutions aperiodic, the projection data were padded with 
zeros to make each projection 256 elements long. 

3.4 Reconstruction from Fan Projections 

The theory in the preceding subsections dealt with reconstructing images 
from their parallel projections such as those shown in Fig. 3. 1. In generating 
these parallel data a source-detector combination has to linearly scan over the 
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Fig. 3.17: (a) Reconstruction 
obtained by using the filter shown 
in Fig. 3 J 6. The 127 rays in the 
projection were zero-padded so 
that each projection was 256 
elements long. The unit sample 
response h(nr) was used with n 
ranging from - 128 to 127, 
yielding 256 points for this 
function. The number of 
projections was 100 and the 
display matrix size is 128 x 128. 
(b) A numerical comparison of 
they= -0.605 line of the 
reconstruction in (a) with the true 
values. Note that the dishing and 
dc shift artifacts visible in Fig. 
3. 13 have disappeared. (From 
[Ros82].) 
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length of a projection, then rotate through a certain angular interval, then scan 
linearly over the length of the next projection, and so on. This usually results 
in times that are as long as a few minutes for collecting all the data. A much 
faster way to generate the line integrals is by using fan beams such as those 
shown in Fig. 3.3. One now uses a point source of radiation that emanates a 
fan-shaped beam. On the other side of the object a bank of detectors is used to 
make all the measurements in one fan simultaneously. The source and the 
entire bank of detectors are rotated to generate the desired number of fan 
projections. As might be expected, one has to pay a price for this simpler and 
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faster method of data collection; as we will see later the simple backprojec- 
tion of parallel beam tomography now becomes a weighted backprojection. 

There are two types of fan projections depending upon whether a 
projection is sampled at equiangular or equispaced intervals. This difference 
is illustrated in Fig. 3.18. In (a) we have shown an equiangular set of rays. If 
the detectors for the measurement of line integrals are arranged on the 
straight line D\D 29 this implies unequal spacing between them. If, however, 
the detectors are arranged on the arc of a circle whose center is at S, they may 
now be positioned with equal spacing along this arc (Fig. 3.18(b)). The 
second type of fan projection is generated when the rays are arranged such 
that the detector spacing on a straight line is now equal (Fig. 3.18(c)). The 
algorithms that reconstruct images from these two different types of fan 
projections are different and will be separately derived in the following 
subsection. 

3.4.1 Equiangular Rays 

Let Rp(y) denote a fan projection as shown in Fig. 3.19. Here 0 is the 
angle that the source S makes with a reference axis, and the angle y gives the 
location of a ray within a fan. Consider the ray SA. If the projection data 
were generated along a set of parallel rays, then the ray SA would belong to a 
parallel projection P 9 (t) for 0 and / given by 

0 = 0 + 7 and / = Z>sin7 (70) 

where D is the distance of the source S from the origin O. The relationships 
in (70) are derived by noting that all the rays in the parallel projection at angle 
9 are perpendicular to the line PQ and that along such a line the distance OB is 
equal to the value of /. Now we know that from parallel projections P e (t) we 
may reconstruct /(x, y) by 

/(jc, y)= V \' m P e (t)h(x cos 0+y sin 0-t)dt dd (71) 

Jo J-/ m 

where t m is the value of t for which P e (t) = 0 with | / 1 > t m in all projections. 
This equation only requires the parallel projections to be collected over 180° . 
However, if one would like to use the projections generated over 360°, this 
equation may be rewritten as 

/(*, y) = \ P [ m P e (t)h(x co$$+y sin 6 - /) dt dd. (72) 

Derivation of the algorithm becomes easier when the point (*, y) (marked C 
in Fig. 3.20) is expressed in polar coordinates (r, 0), that is, 

x=rcos<£ >> = rsin<J>. ( 7 ^) 
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Rays at equiangular 
intervals 




(a) 



Detector 
spac ing 
unequal 



Fig. 3.18: Two different types 
of fan beams are shown here* In 

(a) the angle between rays is 
constant but the detector spacing 
is uneven. If the detectors are 
placed along a circle the spacing 
mil then be equal as shown in 

(b) . As shown in (c) the detectors 
can be arranged with constant 
spacing along a line but then the 
angle between rays is not 
constant. (From {Ros82J.) 




2 Da 
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Angular interval 




(C) 



Fig. 3.18: Continued. The expression in (72) can now be written as 

fir, <£) P P" P$(f)h(r cos (0 - 0) - /) dt dO. (74) 

Using the relationships in (70), the double integration may be expressed in 
terms of y and j3, 

f(r,<l>) = - P tf+7 (Z> sin 7)*(r cos (jS + 7-*) 

2 J -7 J -sin" 1 (r m /i>) 

- sin 7)i? cos 7 rf7 rfjS (75) 

where we have used dt dd - D cos y dy d(3. A few observations about this 
expression are in order. The limits - 7 to 2ir - 7 for j3 cover the entire range 
of 360°. Since all the functions of /S are periodic (with period 2tc) these limits 
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may be replaced by 0 and 2tt, respectively. Sin" 1 (t m /D) is equal to the value 
of 7 for the extreme ray SE in Fig. 3.19. Therefore, the upper and lower 
limits for 7 may be written as y m and -y m , respectively. The expression 
P(3 +y (D sin 7) corresponds to the ray integral along SA in the parallel 
projection data P 9 (t). The identity of this ray integral in the fan projection 
data is simply #0(7). Introducing these changes in (75) we get 

/(r, 4,) = ! ( 2 * P m R 0 (y)h(rcos(P + y-<t>)-Dsiny)Dcosydydl3. 

(76) 

In order to express the reconstruction formula given by (76) in a form that 
can be easily implemented on a computer we will first examine the argument 
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Fig. 3.19: An equiangular fan is 
shown here. Each ray is identified 
by its angle 7 from the central 
ray. (From fRos82J.) 



Fig. 3.20: This figure illustrates 
that L is the distance of the pixel 
at location (x, y)from the source 
S; and y is the angle that the 
source-to-pixel line subtends with 
the central ray. (From [Ros82J.) 



of the function h. The argument may be rewritten as 

r cos (0 + y sin 7 

= r cos (0 - 4>) cos 7 - [r sin (0 - <f>) + D] sin 7. (77) 

Let L be the distance from the source 5 to a point (x, y) [or (r, <f>) in polar 
coordinates] such as C in Fig. 3.20. Clearly, L is a function of three 
variables, r, 4> 9 and j3. Also, let 7 ' be the angle of the ray that passes through 
this point (r, <f>). One can now easily show that 

L cos 7' = Z?+r sin (/?-<£) 

Lsin7'=rcos(j8-^). (78) 

Note that the pixel location (r, <£) and the projection angle 0 completely 
determine both L and 7' : 



and 



L(r y <f> 9 |S) = V[Z) + r sin ( 0 - <t>)] 2 + [r cos ( 0 - <A)] 2 (79) 



, r cos (/?-<£) 

7' =tan-» . (80) 

£> + r sin (j3-<£) 
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Using (78) in (77) we get for the argument of h 

rcos (0 + 7-0)-£>sin y = JL sin (7' -7) ( 81 ) 
and substituting this in (76) we get 

/ (r> 0) J- ( 2t V m RfaML sin (7' -7))^ cos 7 dy dfi. (82) 

2 J 0 J -y m 

We will now express the function h(L sin (7' - 7)) in terms of h{t). Note 
that h{t) is the inverse Fourier transform of | w| in the frequency domain: 

A(0=r \w\e JUwt dw. (83) 

J — OD 

Therefore, 

h(L sin 7)= |w|e> 2 * wLsin ? </w. (84) 
Using the transformation 

W '=^2 (85) 



we can wnte 



(86) 



h(Lsiny)=(— Y — V P \w'\e»'»'y dw' 
\Ls\ny/ J -» 

= (-^-Va(7). < 87 > 

\L sin 7 / 

Therefore, (82) may be written as 

/(r, <t>) = I'' P" ^(T)*(7 ' " T)^ cos 7 tf T 40 (88) 

Jo L l J -y m 



where 



8(y) = U-r-Y *(7). ( 89 > 
2 \ sin 7 / 

For the purpose of computer implementation, (88) may be interpreted as a 
weighted filtered backprojection algorithm. To show this we rewrite (88) as 
follows: 

/(r,*)=f o *^a(T')^ < 90 > 
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where 
and where 

R^(y) = R 0 (y) • D • cos 7. (92) 
This calls for reconstructing an image using the following three steps: 
Step 1: 

Assume that each projection Rp(y) is sampled with sampling interval a. 
The known data then are R 0i (na) where n takes integer values, ft are the 
angles at which projections are taken. The first step is to generate for 
each fan projection R fii (na) the corresponding R'^(na) by 

Rp {na) = Rto{na) * D • cos net. (93) 

Note that n = 0 corresponds to the ray passing through the center of the 
projection. 

Step 2: 

Convolve each modified projection R%(na) with g(na) to generate the 
corresponding filtered projection: 

Qe i {na) = R; i (na)*g(na). (94) 

To perform this discrete convolution using an FFT program the function 
Rfyina) must be padded with a sufficient number of zeros to avoid 
interperiod interference artifacts. The sequence g(not) is given by the 
samples of (89): 

g{na) J (j!lJ) 2 h{na) . (95) 
2 \sin na/ 

If we substitute in this the values of h(not) from (61), we get for the 
discrete impulse response 



r 

1 



g(nct) = < 



w = 0 



8a 2 ' 

0, n is even (96) 



/ a V 

( ) , n is odd. 

\ Tra sin na / 



Although, theoretically, no further filtering of the projection data than 
that called for by (94) is required, in practice superior reconstructions 
are obtained if a certain amount of smoothing is combined with the 
required filtering: 

Q 0i (na) = /?;,(wa)*g(/ia)**(*a) (97) 
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Fig. 3.21: While the filtered 
projections are backprojected 
along parallel lines for the 
parallel beam case (a), for the 
fan beam case the backprojection 
is performed along converging 
lines (bh (c) This figure 
illustrates the implementation step 
that in order to determine the 
backprojected value at pixel (jf, 
y), one must first compute y ' for 
that pixel. (From [Ros82].) 
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s 



(c) 

Fig. 3.21: Continued. where k(na) is the impulse response of the smoothing filter. In the 

frequency domain implementation this smoothing filter may be a simple 
cosine function or a Hamming window. 

Step 3: 

Perform a w/g/tf erf backprojection of each filtered projection along the 
fan. Since the backprojection here is very different from that for the 
parallel case, we will explain it in some detail. For the parallel case the 
filtered projection is backprojected along a set of parallel lines as shown 
in Fig. 3.21(a). For the fan beam case the backprojection is done along 
the fan (Fig. 3.21(b)). This is dictated by the structure of (90): 

where y' is the angle of the fan beam ray that passes through the point 
(jc, y) and A/J = 2tt/M. For ft chosen in Fig. 3 .2 1 (c) in order to find the 
contribution of Q^(y) to the point (x, y) shown there one must first find 
the angle, y\ of the ray SA that passes through that point (x, y). 
Qfyiy') will then be contributed from the filtered projection at ft to the 
point (x, y) under consideration. Of course, the computed value of 7' 
may not correspond to one of noc for which Q^(na) is known. One must 
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then use interpolation. The contribution Q 0i (y f ) at the point (x, y) must 
then be divided by L 2 where L is the distance from the source S to the 
point (x, y). 

This concludes our presentation of the algorithm for reconstructing projection 
data measured with detectors spaced at equiangular increments. 

3.4.2 Equally Spaced Collinear Detectors 



Fig. 3.22: For the case of 
equispaced detectors on a straight 
line, each projection is denoted 
by the function R & (s). (From 
[Ros82]J 



Let R$(s) denote a fan projection as shown in Fig. 3.22, where s is the 
distance along the straight line corresponding to the detector bank. The 
principal difference between the algorithm presented in the preceding 
subsection and the one presented here lies in the way a fan projection is 
represented, which then introduces differences in subsequent mathematical 
manipulations. Before, fan projections were sampled at equiangular intervals 
and we represented them by Rp(y) where y represented the angular location 
of a ray. Now we represent them by Rp(s). 

Although the projections are measured on a line such as D X D 2 in Fig. 3.22, 
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for theoretical purposes it is more efficient to assume the existence of an 
imaginary detector line D[D{ passing through the origin. We now associate 
the ray integral along SB with point A on D[ D{ , as opposed to point B on 
D X D 2 . Thus in Fig. 3.23 we will associate a fan projection R${s) with the 
imaginary detector line D[D{ . Now consider a ray SA in the figure; the 
value of s for this ray is the length of OA. //"parallel projection data were 
generated for the object under consideration, the ray SA would belong to a 
parallel projection P e (t) with 6 and / as shown in the figure. The relationship 
between /3 and t for the parallel case is given by 

t = scosy 0 = j3 + 7 

sD . s 

t= 0^ + tan- 1 - (99) 

V^T^2 D 

where use has been made of the fact that angle AOC is equal to angle OSC, 
and where D is the distance of the source point S from the origin O. 

In terms of the parallel projection data the reconstructed image is given by 
(74) which is repeated here for convenience: 

/( r> = I [ 2T \' m P $ (t)h(r cos (•-*)- /) dt dd (74) 
2 J o J -t m 



s 
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Fig. 3.23: This figure illustrates 
several of the parameters used in 
the derivation of the 
reconstruction algorithm for 
equispaced detectors. (From 
fRos82J.) 



where /(r, </>) is the reconstructed image in polar coordinates. Using the 
relationships in (99) the double integration may be expressed as 



fir, <A)=- , P^A — - 

(100) 

where we have used 

* rf '" (DW)M * rf * (101) 

In (100) 5 m is the largest value of 5 in each projection and corresponds to t m 
for parallel projection data. The limits -tan -1 (s m /D) and 27r-tan~ I (s m / 
D) cover the angular interval of 360°. Since all functions of j3 in (100) are 
periodic with period 2tt, these limits may be replaced by 0 and 2tt, 
respectively. Also, the expression 



(102) 



corresponds to the ray integral along SA in the parallel projection data P${t). 
The identity of this ray integral in the fan projection data is simply Rq(s). 
Introducing these changes in (100) we get 

f(r, <t>)= X - J** j*" Re(s)h (r cos (jS+tan" 1 ^-0^ 



Ds \ 



D i 

ds (103) 



(D 2 +s 2 ) v2 



In order to express this formula in a filtered backprojection form we will first 
examine the argument of h. The argument may be written as 



r cos ^S + tan -1 ^-<£^ - 



Ds 



\fb~ 2 ~+s 1 



D s 
= rcos (0-<£) -(£> + /• sin (/3-0))— = . (104) 

y/D 2 + s 2 y/D 2 + s 2 

We will now introduce two new variables that are easily calculated in a 
computer implementation. The first of these, denoted by U, is for each pixel 
(x, y) the ratio of SP (Fig. 3.24) to the source-to-origin distance. Note that 
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Fig. 3.24: For a pixel at the 
polar coordinates (r, <f>) the 
variable V is the ratio of the 
distance SP, which is the 
projection of the source to pixel 
line on the central ray, to the 
source-to-center distance. 
(Adapted from [Ros82].) 



SP is the projection of the source to pixel distance SE on the central ray. Thus 



t/(r, 4>,]8) = - 



D 

£> + r sin ()8-<£) 
= D 



(105) 



(106) 



The other parameter we want to define is the value of s for the ray that passes 
through the pixel (r, 4>) under consideration. Let s' denote this value of s. 
Since s is measured along the imaginary detector line D[Di , it is given by 
the distance OF. Since 



we have 



s' EP 
SO SP 

r cos (&-<(>) 
Z> + r sin (/3-<£) 



(107) 



(108) 
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Equations (106) and (108) can be utilized to express (104) in terms of U and 

s': 

( s \ Ds s'UD sUD 

rcos fl + tan- 1 --* )-- = =- . • (109) 

\ D J y/D 2 +s 2 y/D 2 +s 2 V£> 2 +s 2 

Substituting (109) in (103), we get 

(110) 

We will now express the convolving kernel h in this equation in a form closer 
to that given by (61). Note that, nominally, h(t) is the inverse Fourier 
transform of | w| in the frequency domain: 

A(0=r \w\e*""dw. (Ill) 

J - as 

Therefore, 

h\(s'-s)-^=] = r \w\e»™«' -'W^ 51 ^) dw. (112) 
Using the transformation 

w' = w UD (113) 
y/D 2 + s 2 

we can rewrite (112) as follows: 

* b'- s) ^\ - w n. i-'i^"-'"' (»«> 
=W* (j, - s) - <|15) 

Substituting this in (110) we get 

/<'• ♦> - r tp r. ^ w *< j ' -*> rm? oi6) 

where 

g(5) = ^/l(5). 017) 

For the purpose of computer implementation, (1 16) may be interpreted as a 
weighted filtered backprojection algorithm. To show this we rewrite (1 16) as 
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follows: 



p VZ) 2 + 5 2 



where 



and 



(118) 



(119) 



(120) 



Equations (118) through (120) suggest the following steps for computer 
implementation: 



Step 1: 



Assume that each projection R & (s) is sampled with a sampling interval of 
a. The known data then are R^(na) where n takes integer values with n 
= 0 corresponding to the central ray passing through the origin; ft are 
the angles for which fan projections are known. The first step is to 
generate for each fan projection R^(na) the corresponding modified 
projection R^(na) given by 



R£ i (na) = Rfi i (na) 



\/D z + n 2 a 2 



(121) 



Step 2: 



Convolve each modified projection R^(na) with g{na) to generate the 
corresponding filtered projection: 

Q 0i (na)^R Pi (na) * g(na) (122) 

where the sequence g(na) is given by the samples of (1 17): 

g(na)=\h{na). (123) 

Substituting in this the values of h(na) given in (61) we get for the 
impulse response of the convolving filter: 
r 

— , n=0 
8a 2 



g(na)=< 



0, 



1 



n even 
n odd. 



(124) 
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When the convolution of (122) is implemented in the frequency domain 
using an FFT algorithm the projection data must be padded with a 
sufficient number of zeros to avoid distortion due to interperiod 
interference. 

In practice superior reconstructions are obtained if a certain amount of 
smoothing is included with the convolution in (122). If k(na) is the 
impulse response of the smoothing filter, we can write 

Q^(na) = R 0i (na) * g(na) * k(na). (125) 

In a frequency domain implementation this smoothing may be achieved 
by a simple multiplicative window such as a Hamming window. 

Step 3: 

Perform a weighted backprojection of each filtered projection along the 
corresponding fan. The sum of all the backprojections is the recon- 
structed image 

to»-"%ln£j3i<W (126) 

where U is computed using (106) and s' identifies the ray that passes 
through (x, y) in the fan for the source located at angle ft. Of course, 
this value of s' may not correspond to one of the values of na at which 
Qp. is known. In that case interpolation is necessary. 



3.4.3 A Re-sorting Algorithm 

We will now describe an algorithm that rapidly re-sorts the fan beam 
projection data into equivalent parallel beam projection data. After re-sorting 
one may use the filtered backprojection algorithm for parallel projection data 
to reconstruct the image. This fast re-sorting algorithm does place constraints 
on the angles at which the fan beam projections must be taken and also on the 
angles at which projection data must be sampled within each fan beam 
projection. 

Referring to Fig. 3.19, the relationships between the independent variables 
of the fan beam projections and parallel projections are given by (70): 

t = D sin y and $ = fi + y. (127) 

If, as before, R$(y) denotes a fan beam projection taken at angle 0, and P e (t) 
a parallel projection taken at angle 0, using (127) we can write 

Re(y) = Pe + y(Dsiny). (128) 

Let A/3 denote the angular increment between successive fan beam 
projections, and let Ay denote the angular interval used for sampling the fan 
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beam projections. We will assume that the following condition is satisfied: 

AjS = A 7 = a. (129) 

Clearly then /? and 7 in (128) are equal to met and not, respectively, for some 
integer values of the indices m and n. We may therefore write (128) as 



(na) = P im+n) a(D sin na). 



(130) 



This equation serves as the basis of a fast re-sorting algorithm. It expresses 
the fact that the /ith ray in the mth radial projection is the /ith ray in the (m + 
ri)th parallel projection. Of course, because of the sin not factor on the right- 
hand side of (130), the parallel projections obtained are not uniformly 
sampled. This can usually be rectified by interpolation. 

3*5 Fan Beam Reconstruction from a Limited Number of Views 

Simple geometrical arguments should convince the reader that parallel 
projections that are 180° apart, P $ (t) and Pe+m<0, are mirror images of 
each other. That is, 

P*(0 = ^ + iso°(-0 031) 

and thus it is only necessary to measure the projections of an object for angles 
from 0 to 180°. 

We can extend this result by noting that an object is completely specified if 
the ray integrals of the object are known for 

0 o <0<0o+18O° (132) 

and 

-Ux^'max U33) 



Fig. 3.25: As shown in this 
figure, each line integral can be 
thought of as a single point in the 
Radon transform of the object. 
Each line integral is identified by 
its distance from the origin and 
its angle. 



where is large enough so that each projection is at least as wide as the 
object at its widest. If each ray integral is represented as a point in a polar 
coordinate system (t, 6) as shown in Fig. 3.25 then a complete set of ray 
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(a) (b) (c) 

integrals will completely fill a disk of radius . This is commonly known as 
the Radon transform or a sinogram and is shown for the Shepp and Logan 
phantom both in polar and rectangular coordinates in Fig. 3.26. 

These ideas can also be extended to the fan beam case. From Fig. 3.27 we 
see that two ray integrals represented by the fan beam angles 71) and (j3 2> 
72) are identical provided 

jSi-7i = &-72+180 o (134) 

and 

7i= "72. 035) 

With fan beam data the coordinate transformation 

i = D sin 7 

0 = /3 + 7 O 36 ) 

maps the (0, 7) description of a ray in a fan into its Radon transform 
equivalent. This transformation can then be used to construct Fig. 3.27, 
which shows the data available in Radon domain as the projection angle 0 
varies between 0 and 180° with a fan angle of 40° (7^ = 20°). 

Recall that points in Radon space that are periodic with respect to the 
intervals shown in (132) and (133) represent the same ray integral. Thus the 
data in Fig. 3.28 for angles 0 > 180° and t > 0 are equal to the Radon data 
for $ < 0 and / < 0. These two regions are labeled A in Fig. 3.28. On the 
other hand, the regions marked B in Fig. 3.28 are areas in the Radon space 
where there are no measurements of the object. To cover these areas it is 
necessary to measure projections over an additional 27 m degrees as shown in 



Fig. 3.26: An object and its 
Radon transform are shown here. 
The object in (a) is used to 
illustrate the short scan algorithm 
developed by Parker fPar82aJ, 

(b) shows the Radon transform 
in rectangular coordinates, while 

(c) represents the Radon 
transform in polar coordinates. 
(Reprinted with permission from 
[Par82a], [Par82b].) 
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Fig, 3.27: Rays in two fan 
beams will represent the same line 
integral if they satisfy the 
relationship fi t - y t = & - 72 




Fig. 3.28: Collecting projections 
over J80° gives estimates of the 
Radon transform between the 
curved lines as shown on the left. 
The curved lines represent the 
most extreme projections for a 
fan angle ofy m . On the right is 
shown the available data in the 
0-y coordinate system used in 
describing fan beams. In both 
cases the region marked A 
represents the part of the Radon 
transform where two estimates 
are available. On the other hand, 
for 180* of projections there are 
no estimates of the Radon 
transform in the regions marked 
B. 
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Fig. 3.29: If projections are 
gathered over an angle of 180° + 
2y m then the data illustrated are 
available. Again on the left is 
shown the Radon transform while 
the right shows the available data 
in the ($-y coordinate system. The 
line integrals in the shaded 
regions represent duplicate data 
and these points must be 
gradually weighted to obtain 
good reconstructions. 



Fig. 3.29. Thus it should be possible to reconstruct an object using fan beam 
projections collected over 180 + 2y m degrees. 

Fig. 3.30 shows a "perfect" reconstruction of a phantom used by Parker 
[Par82a], [Par82b] to illustrate his algorithm for short scan or "180 degree 
plus" reconstructions. Projection data measured over a full 360° of )3 were 
used to generate the reconstruction. 

It is more natural to discuss the projection data overlap in the (/J, 7) 
coordinate system. We derive the overlap region in this space by using the 
relations in (134) and (135) and the limits 

0<^<180°+2 7w 

0<j3 2 <180° + 2 7m . (137) 

Substituting (134) and (135) into the first equation above we find 

0 < 0 2 - 2 72 + 1 80° < 1 80 0 + 2y m (138) 

and then by rearranging 

- 1 80 D + 2 72 < ft * 2y m -2y 2 . (139) 

Substituting the same two equations into the second inequality in (137) we 
find 

0< 13, - 27, - 180° < 180° 4- 2y m (140) 
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Fig. 3.30: This figure shows a 
reconstruction using 360* of fan 
beam projections and a standard 
filtered backprojection algorithm, 
(Reprinted with permission from 
[Par82a] t [Par82b].) 



and then by rearranging 

180° + 2 7l < ft < 360° + 2y m + 2y } . 



(141) 



Since the fan beam angle, 7, is always less then 90°, the overlapping regions 
are given by 



O<0 2 <2 7m + 272 



(142) 



Fig. 3.31: This reconstruction 
was generated with a standard 
filtered backprojection algorithm 
using 220° of projections. The 
targe artifacts are due to the lack 
of data in some regions of the 
Radon transform and duplicate 
data in others. (Reprinted with 
permission from [Par82a], 
[Par82b].) 



and 



180° + 27i<j8 1 <180° + 27„ 



(143) 



as is shown in Fig. 3.29. 

If projections are gathered over an angle of 180° + 27 m and a 
reconstruction is generated using the standard fan beam reconstruction 
algorithms described in Section 3.4, then the image in Fig. 3.31 is obtained. 
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In this case a fan angle of 40° (7^ = 20) was used. As described above, the 
severe artifacts in this reconstruction are caused by the double coverage of the 
points in region B of Fig. 3.28. 

One might think that the reconstruction can be improved by setting the data 
to zero in one of the regions of overlap. This can be implemented by 
multiplying a projection at angle )8, p 0 (y) 9 by a one-zero window, Wpiy), 
given by 

f0 0</3<27m + 2 7 (144) 
\^ 1 elsewhere. 

As shown by Naparstek [Nap80] using this type of window gives only a small 
improvement since streaks obscure the resulting image. 

While the above filter function properly handles the extra data, better 
reconstructions can be obtained using a window described in [Par82a]. The 
sharp cutoff of the one-zero window adds a large high frequency component 
to each projection which is then enhanced by the | w | filter that is used to filter 
each projection. 

More accurate reconstructions are obtained if a "smoother" window is 
used to filter the data. Mathematically, a "smooth'* window is both 
continuous and has a continuous derivative. Formally, the window, ^(7), 
must satisfy the following condition: 

^,(7i) + w ft (72)=l < 145 ) 

for (j3i, 71) and (j3 2 > 72) satisfying the relations in (134) and (135), and 

w 0 (7) = 0 (146) 

and 

MW +27/n = 0. (147) 

To keep the filter function continuous and * Smooth' ' at the boundary between 
the single and double overlap regions the following constraints are imposed 
on the derivative of ^(7): 



3^(7) 



and 

3wg(7) 



= 0 048) 

0 = 2 7m + 2 7 



= 0. (149) 

0=180° +2? 
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Fig. 3J2: Using a weighting 
function that minimizes the 
discontinuities in the projection 
this reconstruction is obtained 
using 220* of projection data. 
(Reprinted with permission from 
[Par82a], [Par82b].) 



One such window that satisfies all of these conditions is 



sin' 
1, 



r 4^-1 

LTm-TJ 



0</3<2 7m -2 7 



sin 2 [ 45 o 180°-f2 7m /? 1 ^ 180 °-2 7 <i8<180' + 2 7/ , 
L Y + 7m J 



(150) 

A reconstruction using this weighting function is shown in Fig. 3.32. From 
this image we see that it is possible to eliminate the overlap without 
introducing errors by using a smooth window. 



3.6 Three-Dimensional Reconstructions 1 



The conventional way to image a three-dimensional object is to illuminate 
the object with a narrow beam of x-rays and use a two-dimensional 
reconstruction algorithm. A three-dimensional reconstruction can then be 
formed by illuminating successive planes within the object and stacking the 
resulting reconstructions. This is shown in Fig. 3.33. 

A more efficient approach, to be considered in this section, is a 
generalization of the two-dimensional fan beam algorithms presented in 
Section 3.4.2. Now, instead of illuminating a slice of the object with a fan of 
x-rays, the entire object is illuminated with a point source and the x-ray flux is 
measured on a plane. This is called a cone beam reconstruction because the 



1 We are grateful for the help of Barry Roberts in the preparation of this material. 
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object 



reconstructed slices 
of the object 



Fig. 3,33: A three-dimensional 
reconstruction can be done by 
repetitively using two-dimensional 
reconstruction algorithms at 
different heights along the z-axis. 
(From !Kak86].) 



rays form a cone as illustrated in Fig. 3.34. Cone beam algorithms have been 
studied for use with Mayo Clinic's Digital Spatial Reconstructor (DSR) 
[Rob83] and Imatron's Fifth Generation Scanner [Boy83J. 

The main advantage of cone beam algorithms is the reduction in data 
collection time. With a single source, ray integrals are measured through 
every point in the object in the time it takes to measure a single slice in a 
conventional two-dimensional scanner. The projection data, R$(t, r), are 
now a function of the source angle, 0, and horizontal and vertical positions on 
the detector plane, / and r. 



3.6.1 Three-Dimensional Projections 



A ray in a three-dimensional projection is described by the intersection of 
two planes 

/=jccos0+j>sin0 (151) 



r= -(-x sin $+y cos 6) sin 7 + z cos 7. 



(152) 



A new coordinate system (/, 5, r) is obtained by two rotations of the (x, y, z)- 
axis as shown in Fig. 3 .35. The first rotation, as in the two-dimensional case, 
is by 6 degrees around the z-axis to give the (/, s, z)-axes. Then a second 
rotation is done out of the (/, .s)-plane around the f-axis by an angle of 7. In 
matrix form the required rotations are given by 



1 0 

0 cos 7 
0 - sin 7 



0 
sin 7 
cos 7 



cos 6 sin 9 
- sin 0 cos 6 
0 0 



0 
0 
1 



x 

y 
z 



(153) 



A three-dimensional parallel projection of the object /is expressed by the 
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/ Detection 
/ plane 



Fig. 3,34: In cone beam 
projections the detector measures 
the x-ray flux over a plane. By 
rotating the source and detector 
plane completely around the 
object all the data necessary for a 
three-dimensional reconstruction 
can be gathered in the time a 
conventional fan beam system 
collects the data for its 
two-dimensional reconstruction. 
(From fKak86J.) 



following integral: 



(154) 



Note that four variables are being used to specify the desired ray; (f, 0) 
specify the distance and angle in the x-y plane and (r, 7) in the s-z plane. 

In a cone beam system the source is rotated by (3 and ray integrals are 
measured on the detector plane as described by /?/?(/>', D- To fi nd the 
equivalent parallel projection ray first define 



p'D. 



'so 



DE 



rAo 



(155) 



DE 



as was done in Section 3.4.2. Here we have used £>so to indicate the distance 
from the center of rotation to the source and A>e to indicate the distance from 
the center of rotation to the detector. For a given cone beam ray, R$(p, f), 
the parallel projection ray is given by 



so 



0 = jS + tan-' (p/Dso) 
where / and 0 locate a ray in a given tilted fan, and similarly 



(156) 



(157) 



(158) 
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Fig. 3.35: To simplify the 
discussion of the cone beam 
reconstruction the coordinate 
system is rotated by the angle of 
the source to give the (s, t)-axis. 
The r-axis is not shown but is 
perpendicular to the t- and s-axes. 
(From [Kak86D 



T = tan-'(f/^$o). 059) 

were r and 7 specify the location of the tilted fan itself. 

The reconstructions shown in this section will use a three-dimensional 
version of the Shepp and Logan head phantom. The two-dimensional ellipses 
of Table 3.1 have been made ellipsoids and repositioned within an imaginary 
skull. Table 3,2 shows the position and size of each ellipse and Fig. 3.36 
illustrates their position. 

Because of the linearity of the Radon transform, a projection of an object 
consisting of ellipsoids is just the sum of the projection of each individual 



Table 3.2: Summary of parameters for three-dimensional tomography simulations. 





Coordinates of 




Rotation 


Gray 




the Center 


Axis Lengths 


Angle 0 


Level 


Ellipsoid 


{x, y, Z) 


(A, B, O 


(deg) 


P 


a 


(0, 0, 0) 


(0.69, 0.92, 0.9) 


0 


2.0 


b 


(0, 0, 0) 


(0.6624, 0.874, 0.88) 


0 


-0.98 


c 


(-0.22, 0, -0.25) 


(0.41, 0.16, 0.21) 


108 


-0.02 


d 


(0.22, 0, -0.25) 


(0.31, 0.11, 0.22) 


72 


-0.02 


e 


(0,0.1, -0.25) 


(0.046, 0.046, 0.046) 


0 


0.02 


f 


(0, 0.1, -0.25) 


(0.046, 0.046, 0.046) 


0 


0.02 


g 


(-0.8,-0.65, -0.25) 


(0.046, 0.023, 0.02) 


0 


0.01 


h 


(0.06, -0.065, -0.25) 


(0.046, 0.023, 0.02) 


90 


0.01 




(0.06, -0.105, 0.625) 


(0.56, 0.04, 0.1) 


90 


0.02 


j 


(0, 0.1, -0.625) 


(0.056, 0.056, 0.1) 


0 


-0.02 
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(a) 



y 




y 



Fig. 3.36: A three-dimensional 
version of the Shepp and Logan 
head phantom is used to test the 
cone beam reconstruction 
algorithms in this section, (a) A 
vertical slice through the object 
illustrating the position of the 
two reconstructed planes, (b) An 
image at plane B (z = -0.25) 
and (c) an illustration of the level 
of each of the ellipses, (d) An 
image at plane A (z = 0.625) and 
(e) an illustration of the gray 
levels with the slice. (From 
[Kak86].) 




(d) 




-►x 



(e) 
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ellipsoid. If the ellipsoid is constant and described by 



( x 2 y 2 z 2 

ABC (16Q) 



0 otherwise 
then its projection on the detector plane is written 

P $ y (/, r) = 2(> f BC \ a 2 ( 0 9 7 ) _ /2 (C 2 cos 2 + (fl 2 cos 2 Q + A 2 sin 2 fl) §in 2 ?) 

a 2 (0,y) i 

/7 + cos(47)\ 
-r 2 (y4 2 cos 2 0 + B 2 sin 2 0)f 

]l/2 
(161) 

where 

a 2(0, 7 )^c 2 (5 2 sin 2 6 + A 2 cos 2 0) cos 2 y + A 2 B 2 sin 2 7. (162) 
If the tilt angle, 7, is zero then (161) simplifies to (5). 

3.6.2 Three-Dimensional Filtered Backprojection 

We will present a filtered backprojection algorithm based on analyses 
presented in [Fel84] and [Kak86]. The reconstruction is based on filtering and 
backprojecting a single plane within the cone. In other words, each elevation 
in the cone (described by z or f ) is considered separately and the final three- 
dimensional reconstruction is obtained by summing the contribution to the 
object from all the tilted fan beams. 

The cone beam algorithm sketched above is best derived by starting with 
the filtered backprojection algorithm for equispatial rays. In a three- 
dimensional reconstruction each fan is angled out of the source-detector 
plane of rotation. This leads to a change of variables in the backprojection 
algorithm. 

First consider the two-dimensional fan beam reconstruction formula for the 
point (r, <}>): 

Sir. <f>) = \ P i P *fi(P)*(P' -P) -?== d P d $ < 163 > 
2 Jo U 1 J-* Vz^+z? 2 

, Aso'cos (P-<f>) f / v 

p = 



Aso + rsin (<S-0) 



*(/>)=[ |«|e'"-*</co (164) 



Dso + r sin (0-<f>) 

U(r, <t>, J3) = -^ „ . (165) 

^so 
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Equation (163) is the same as (116), except that we have now used different 
names for some of the variables. To further simplify this expression we will 
replace the (r, <f>) coordinate system by the rotated coordinates (/, s). Recall 
that (t, s) is the location of a point rotated by the angular displacement of the 
source-detector array. The expressions 

f=jtcos p+y sin 0 s= -jc sin j8 + y cos 0 (166) 

x=r cos 0 y = r sin <t>, (167) 

lead to 

p . m D*L U{x>y ^ )= ^Zl. (168) 
The fan beam reconstruction algorithm is now written as 

(169) 

In a cone beam reconstruction it is necessary to tilt the fan out of the plane 
of rotation; thus the size of the fan and the coordinate system of the 
reconstructed point change. As shown in Fig. 3.37 a new coordinate system 
(f, s) is defined that represents the location of the reconstructed point with 
respect to the tilted fan. Because of the changing fan size both the source 
distance, D SQ , and the angular differential, j8, change. The new source 
distance is given by 



Z 




s 



Fig. 3.37: The (£ s) coordinate 
system represents a point in the 
object with respect to a tilted fan 
beam. (From [Kak86U 
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where f is the height of the fan above the center of the plane of rotation. In 
addition, the increment of angular rotation becomes 



Dndfi^D^dfi' dp' =T==- (171) 

so^" ^ 

Substituting these new variables, D^t for Dso and <//3' for <//3, and writing 
the projection data as Rp>(p, f), (169) becomes 

■ f fy'(p> m (i^.-p) 7#=i * < 172 > 

J — V^so"* /VZ>so+P 2 

To return the reconstruction to the original (/, s, z) coordinate system we 
substitute 

s s f Z 

t=t, = , = (173) 

Dso Dso Dso-s 



and (170) and (171) to find 

1 f2» 



1 DIq 



(Dso-s) 2 

• J' *<a o* 7== * dfi ' (174) 

The cone beam reconstruction algorithm can be broken into the following 
three steps: 

Step 1: 

Multiply the proje ction data, Rpip, f), by the function (Dso/ 
^so + J* + P 2 ) to find R fo> f> : 

*;(p.0--=%== D. 075) 

Step 2: 

Convolve the weighted projection /^(/>, f) with A(/?)/2 by multiplying 
their Fourier transforms with respect to /?. Note this convolution is done 
independently for each elevation, f. The result, Qe(p, f), is written 

QeiP, t) = ^(P,f)*5*(P). ( 176 > 

Step 3: 

Finally, each weighted projection is backprojected over the three- 
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dimensional reconstruction grid: 

The two arguments of the weighted projection, Qp, represent the 
transformation of a point in the object into the coordinate system of the 
tilted fan shown in Fig. 3.37. 

Only those points of the object that are illuminated from all directions can 
be properly reconstructed. In a cone beam system this region is a sphere of 
radius D so sin (T m ) where T m is half the beamwidth angle of the cone. 
Outside this region a point will not be included in some of the projections and 
thus will not be correctly reconstructed. 

Figs. 3.38 and 3.39 show reconstructions at two different levels of the 
object described in Fig. 3.36. In each case 100 projections of 127 x 127 
elements were simulated and both a gray scale image of the entire plane and a 
line plot are shown. The reconstructed planes were at z = 0.625 and z = 
-0.25 planes and are marked as Plane A and Plane B in Fig. 3.36. 

In agreement with [Smi85], the quality of the reconstruction varies with the 
elevation of the plane. On the plane of rotation (z - 0) the cone beam 
algorithm is identical to a equispatial fan beam algorithm and thus the results 
shown in Fig. 3.38 are quite good. Farther from the central plane each point 
in the reconstruction is irradiated from all directions but now at an oblique 
angle. As shown in Fig. 3.39 there is a noticeable degradation in the 
reconstruction. 

3.7 Bibliographic Notes 

The current excitement in tomographic imaging originated with Houns- 
field's invention [Hou72] of the computed tomography (CT) scanner in 1972, 
which was indeed a major breakthrough. His invention showed that it is 
possible to get high-quality cross-sectional images with an accuracy now 
reaching one part in a thousand in spite of the fact that the projection data do 
not strictly satisfy theoretical models underlying the efficiently implement- 
able reconstruction algorithms. (In x-ray tomography, the mismatch with the 
assumed theoretical models is caused primarily by the polychromaticity of the 
radiation used. This will be discussed in Chapter 4.) His invention also 
showed that it is possible to process a very large number of measurements 
(now approaching a million) with fairly complex mathematical operations, 
and still get an image that is incredibly accurate. The success of x-ray CT has 
naturally led to research aimed at extending this mode of image formation to 
ultrasound and microwave sources. 

The idea of filtered backprojection was first advanced by Bracewell and 
Riddle [Bra67] and later independently by Ramachandran and Lakshminaray- 
anan [Ram71]. The superiority of the filtered backprojection algorithm over 
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(a) 



Fig. 3.38: (a) Cone beam 
algorithm reconstruction of plane 
B in Fig. 3.36. (b) Plot of the y 
= - 0.605 line in the 
reconstruction compared to the 
original. (From [Kak86].) 
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(3) 



Fig- 3.39: (a) Cone beam 
algorithm reconstruction of plane 
A in Fig. 3.36. (b) Plot of the y 
= -0.105 line in the 
reconstruction compared to the 
original. (From fKak86J.) 
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the algebraic techniques was first demonstrated by Shepp and Logan [She74]. 
Its development for fan beam data was first made by Lakshminarayanan 
[Lak75] for the equispaced collinear detectors case and later extended by 
Herman and Naparstek [Her77] for the case of equiangular rays. The fan 
beam algorithm derivation presented here was first developed by Scudder 
[Scu78]. Many authors [Bab77], [Ken79], [Kwo77], [Lew79], [Tan75] have 
proposed variations on the filter functions of the filtered backprojection 
algorithms discussed in this chapter. The reader is referred particularly to 
[Ken79], [Lew79] for ways to speed up the filtering of the projection data by 
using binary approximations and/or inserting zeros in the unit sample 
response of the filter function. Images may also be reconstructed from fan 
beam data by first sorting them into parallel projection data. Fast algorithms 
for ray sorting of fan beam data have been developed by Wang [Wan77], 
Dreike and Boyd [Dre77], Peters and Lewitt [Pet77], and Dines and Kak 
[Din76]. The reader is referred to [Nah81] for a filtered backprojection 
algorithm for reconstructions from data generated by using very narrow angle 
fan beams that rotate and traverse continuously around the object. The 
reader is also referred to [Hor78], [Hor79] for algorithms for nonuniformly 
sampled projection data, and to [Bra67], [Lew78], [Opp75], [Sat80], 
[Tam81] for reconstructions from incomplete and limited projections. Full 
three-dimensional reconstructions have been discussed in [Chi79], [Chi80], 
[Smi85]. We have also not discussed the circular harmonic transform method 
of image reconstruction as proposed by Hansen [Han81a], [Han81b]. 

Tomographic imaging may also be accomplished, although less accurately, 
by direct Fourier inversion, instead of the filtered backprojection method 
presented in this chapter. This was first shown by Brace well [Bra56] for 
radioastronomy, and later independently by DeRosier and Klug [DeR68] in 
electron microscopy and Rowley [Row69] in optical holography. Several 
workers who applied this method to radiography include Tretiak et al. 
[Tre69], Bates and Peters [Bat71], and Mersereau and Oppenheim [Mer74]. 
In order to utilize two-dimensional FFT algorithms for image formation, the 
direct Fourier approach calls for frequency domain interpolation from a polar 
grid to a rectangular grid. For some recent methods to minimize the resulting 
interpolation error, the reader is referred to [Sta81]. Recently Wernecke and 
D'Addario [Wer77] have proposed a maximum-entropy approach to direct 
Fourier inversion. Their procedure is especially applicable if for some reason 
the projection data are insufficient. 
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